arXiv:1506.07360vl [cond-mat.mes-hall] 24 Jun 2015 
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The finite-temperature transport properties of the spinless interacting fermion model coupled to 
non-interacting leads are investigated. Employing the unrestricted time-dependent Hartree-Fock 
(HF) approximation, the transmission probability and the non-linear I-V characteristics are cal¬ 
culated, and compared with available analytical results and with numerical data obtained from a 
Hubbard-Stratonovich decoupling of the interaction. In the weak interaction regime, the HF ap¬ 
proximation reproduces the gross features of the exact I-V characteristics but fails to account for 
subtle properties like the particular power law for the reflected current in the interacting resonant 
level model. 


INTRODUCTION 


Out-of-equilibrium quantum systems have received 
much attention in the past few decades, both experimen¬ 
tally and theoretically [T]. A major goal is to understand 
the transport of charge and energy through molecular or 
nanoscale systems coupled to reservoirs such that volt¬ 
ages and temperature gradients can be applied. While 
there exist powerful numerical and analytical methods 
to calculate ground state and finite temperature proper¬ 
ties of isolated interacting quantum systems, the situa¬ 
tion becomes much more involved when these systems are 
coupled to reservoirs and driven out of equilibrium, even 
in the case when a stationary state is reached [5] . Due to 
these difficulties most of the previous studies have been 
restricted to single site models like the spinless interact¬ 
ing resonant level model (IRLM) or the single impurity 
Anderson model, and the main focus was on the zero tem¬ 
perature I-V characteristics, in particular in the linear 
regime. A variety of methods, both numerical and ana¬ 
lytical, have been applied like, e.g., the time-dependent 
density matrix renormalization group (tdDMRG) [SHI], 
non-equilibrium Green’s functions ismi, the time evolv¬ 
ing block decimation method mm, iterated summation 
of the path integral cni, and time-dependent density 
functional theory (tdDFT) [TTI - tTS] . While in the Meir- 
Wingreen approach [6] the problem is formally solved, for 
interacting models it is generally not possible to evaluate 
the various Green’s functions needed as input without 
further approximations. In the purely numerical meth¬ 
ods there exist severe limitations with respect to size and 
dimensionality of the systems that can be studied, and 
even for single site models the approaches are computa¬ 
tionally very expensive. 

In the present study we utilize the Hartree-Fock ap¬ 
proximation in order to calculate the I-V characteris¬ 
tics of a spinless fermion model that can be seen as a 
generalization of the IRLM to several sites. The ob¬ 
vious advantages of this method are the relatively low 
computational costs compared to numerically exact ap¬ 
proaches, and the great flexibility with respect to dimen¬ 


sionality, system size and geometry, finite temperature, 
and type and range of interactions. On the other hand, it 
is well known that HF calculations for isolated systems 
in the ground state or in thermal equilibrium tend to 
overestimate the appearance of spurious broken symme¬ 
try phases. Therefore it is most important to benchmark 
the results against exact solutions. The purpose of this 
study is to assess the reliability of the HF approach in 
the non-equilibrium setting in comparison with available 
exact results. 


MODEL AND METHODS 


Spinless fermion model 


We consider a one-dimensional model of spinless 
fermions, where Nq central sites (the molecule) with 
nearest-neighbor interaction U are coupled to a left and a 
right lead of AA and N-n non-interacting sites. The hop¬ 
ping parameter Iq is, for simplicity, chosen to be the same 
in the leads and within the molecule, while the coupling 
between the leads and the molecule is described by a hop¬ 
ping parameter P and a nearest neighbor interaction U'. 
The Hamiltonian reads 


H = Hh + HhC + He + Hcr + Hr 


( 1 ) 
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FIG. 1. Illustration of the model of Eq. Q for A^l = A^r = 3 
and Nc = 2. 


where hi = cjc;. The sites a + I and 6 are the first and 
the last sites of the molecule, respectively (see Fig. 1), 
and the system size is A^ = A^l + Nc + A^r. The charge 
currents between the leads and the molecule can be ob¬ 
tained from the equation of motion in the Heisenberg 
picture. Assuming that each fermion carries the charge 
e, the charge current operator from the left lead into the 
molecule is given by 


smaller than the band width of the leads [TH]. However, 
for larger voltages there are significant differences in the 
/-+ characteristics. While in the partitioned scheme the 
current approaches a finite value in the limit H —>■ oo, 
in the second scheme it goes to zero when the potential 
difference exceeds the band width. 

Since it is not clear how to realize a well-defined tem¬ 
perature bias in the second setting, and since we want 
to retain the possibility to treat both voltage and tem¬ 
perature gradients we use the first scheme for our cal¬ 
culations. The merits and shortcomings of both quench 
schemes have been extensively discussed in the literature 

[smiiiin]. 


Non-interacting system 


lL = -e 


dNi^ 

dt 


= -ie[H,Nc\ = -iet'{clc^+i - 6 ^+^), 


(3) 

where A^l is the number of particles in the left lead, and 
h is set to one. Correspondingly, the operator of the 
current from the molecule into the right lead reads 


/r = e-^ = ie[H, iVR] = -iet'- 4+i4)- (4) 

Due to particle number conservation, in a stationary 
state the expectation values of left and right currents 
are identical, (/l) = ( 4 r). 

In order to probe the transport properties of a sys¬ 
tem we drive it out of equilibrium by applying a voltage 
and/or a temperature gradient. There are two different 
schemes that have been proposed in the literature USE] 
to simulate a non-equilibrium condition. In the first set¬ 
ting the molecule is initially considered non-interacting 
and decoupled from the leads, and each of the three iso¬ 
lated subsystems is assumed to be in grand canonical 
thermal equilibrium at its own temperature and chemical 
potential. At time t = 0 the interaction is switched on, 
and the molecule and the leads are connected by adding 
Hcc and i?CR- The time evolution of the system is then 
governed by the full Hamiltonian. For sufficiently long 
leads one expects that after a transient phase a station¬ 
ary state arises where the current is time-independent, 
and the goal is to calculate this steady state current as 
function of the applied temperature and voltage bias. 

In the second scheme one considers a situation where 
initially the leads and the molecule are coupled and in 
thermal equilibrium at the same temperature and chem¬ 
ical potential. At time t = 0 a voltage V is applied by 
adding different potentials ±eV/2 to the leads. Subse¬ 
quently, the time evolution is governed by the full Hamil¬ 
tonian including the voltage bias. The stationary cur¬ 
rents thus obtained coincide with the ones of the first 
quench scheme as long as the applied potential is much 


Generally, the state of a quantum mechanical system 
is described by the density matrix p. In our case each 
subsystem a = L,C,R is initially in thermal equilibrium 
with inverse temperature /3 q, and chemical potential pa, 
and the density matrix reads 

/3(0) = g—/3 c(Ac—M eAc)g—/3 r(Ar—/jr+r) 

(5) ’ 

where Z is determined from the condition Tr{/5(0)} = 1. 
The time-evolution of the density matrix is determined 
through 


m = (6) 

In order to calculate time-dependent expectation values 
of observables, like the current, one has to compute the 
equal-time Green’s function defined as 

Gim = {clnCi) = Tr{p{t)cl,ci}. (7) 

For non-interacting systems, i.e., when the Hamiltonian 

H contains only bilinear combinations of the Fermi op¬ 
erators, G can be obtained from the equation of motion 

|G(<) = -z[i 7 ,G], ( 8 ) 

where the matrix H is defined by 

H = Himc\cm- ( 9 ) 

l,m 

Equation (|^ is valid also for time-dependent but non¬ 
interacting Hamiltonians. Expectation values of arbi¬ 
trary products of Fermi operators can be expressed in 
terms of the Green’s function Gim by applying Wick’s 
theorem. Therefore, once the matrix H is numeri¬ 
cally diagonalized, it is straightforward to calculate time- 
dependent charge currents for arbitrary times. 
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Hartree-Fock approach 

When interactions are taken into account the task of 
determining the time-evolution of the density matrix be¬ 
comes much more involved. Exact solutions for interact¬ 
ing systems out of equilibrium are extremely rare and 
limited to very special points in the parameter space of 
the considered model, e.g., the self dual point U' = 2to of 
the IRLM [20] • A rather simple and physically intuitive 
approximation arises from the Hartree-Fock decoupling 
of the interaction terms in the Hamiltonian, 

hifii+i -^{ni)ni+i + {ni+i)ni - {hi+i){hi) - 

- (cj'+icjcj'q+i -b (c|-HiC;)(cj'c/-Hi)- (10) 

As a result, the time-evolution of the density matrix 
is governed by a non-interacting but time-dependent 
Hamiltonian i?HF, and the equation of motion for the 
one-particle density matrix reads 

—GHF(i) = —*[^^HF(i), Ghf(0]- (11) 

This equation of motion can be solved numerically with 
arbitrary precision using the short-time propagation 

Gupit + At) ~ (12) 


and k = {2 tt/Ni^)v with 0 < v < Ni^/2. Here we have 
assumed for simplicity that the wave functions in both 
leads are continued periodically up to infinity, and that 
IVl = Ar. Numerically it is straightforward to connect 
the plane wave ansatz (14) in the left and right lead by 
solving the Schrddinger equation for the HE Hamiltonian 
within the molecule. Since Ahf itself depends on the 
equal-time Green’s function, Eq. (13) has to be solved 
self-consistently. Numerically this is achieved through it¬ 
eration. Starting with some reasonable guess for 
calculates Ahf using Eq. (10) and from there again G 
via Eq. (13). The procedure is iterated until convergence 
is reached. From the self-consistent solution the steady 
state current can then be calculated using Eq. 

The expression for the stationary current can also be 
cast into the form of the usual Landauer formula 


one 

HF 


1=^1 de{Me)-M{e))T{e). (15) 

where T{ek) = is the transmission probability. In 
contrast to the non-interacting case, here T{e) is not only 
a function of the energy, but depends also on the volt¬ 
age and the temperature of the leads due to the self- 
consistency condition. 


Discrete Hubbard-Stratonovich transformation 


and choosing sufficiently small time steps At. 

So far, we have discussed the time-evolution of the sys¬ 
tem starting from a given initial state in order to ex¬ 
tract the transport properties from the stationary state 
that emerges in the course of time. As an alternative, 
one can use an approach that allows one to calculate 
the stationary currents directly without considering ex¬ 
plicitly the time-evolution. Formally this is achieved by 
the Meir-Wingreen approach mzi, based on the non¬ 
equilibrium Green’s function formalism. For an effec¬ 
tively non-interacting system like Arf their result for 
the I-V relation is equivalent to the Landauer formula 
[HiEg. It has been shown [23] that the stationary state 
Green’s function G^^ can be expressed in terms of scat¬ 
tering states as 


If the number of interacting sites is small, ideally if 
there is only interaction across a single link of neighbor¬ 
ing sites as it is the case for Nc = 2 and U' = 0, it is 
possible to calculate the time-dependent density matrix 
p(t) exactly with moderate numerical effort. The start¬ 
ing point is to write the time evolution operator [/ as a 
product of M short-time evolution operators, 

U{t) = (16) 

with t = MAt. For small At one may further split each 
exponential by dividing H = T + V into non-interacting 
contributions T and interaction terms V. Using the sym¬ 
metric Trotter breakup, one obtains 


'^Ui<^k){l\ka){ka\m), (13) 

a—L,R k 


where \ka) is the eigenstate of Ahf that corresponds to 
an incoming wave from lead a = L,R with wave number 
k, and fai^k) = -|- 1)“^ with Ck = —2tocosk 

is the Fermi function that accounts for the thermal oc¬ 
cupation of the states within each lead. Explicitly, up to 
a normalization constant the plane wave scattering state 
originating in the left lead is given by 


{m\kL) 


^-ikm fgj. ^ g L 

tk e*'='" for m e R ’ 


(14) 


-iHAt _ -i{f/2+V+fj2)At ^ -if At/2 -iVAt -if At/2 

- t/ C fcl 

(17) 

with an error of 0{At^). Finally, each exponential con¬ 
taining V can be replaced by a sum over an Ising variable 
s = ±1 using Hirsch’s discrete Hubbard-Stratonovich de¬ 
coupling [21] 

g—l/2)(rii + i —l/2)At _ iC/At/4 \ ' ^—ias(ni—ni.^i)At 

2 ^ 

s=±l 

(18) 

with aAt = arccos ^ Applying this transfor¬ 

mation to both exponentials in Eq. § one may express 
p{t) as the summation over the 4^ configurations of 2M 
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FIG. 2. Current voltage characteristics of the IRLM for 
U' = 2to and t'/to = 0.3, 0.5. Symbols: HF approximation, 
full curves: exact results from Eq. (19l. The inset shows (for 
t' — 0.5to) the relative deviation A = (Jq — I)/Io from perfect 
conduction, 7o = GoV, Go = e^/h, on a double logarithmic 
scale. Symbols: HF data, solid line: A oc 


FIG. 4. HF transmission coefficient of the IRLM at zero volt¬ 
age t' — O.lSto and U'/to = 0, 0.5,1,1.5 (from bottom to top). 
The inset shows the broadening F/Fo of the peak as function 
of U' compared to the broadening of the exact spectral func¬ 
tion of the zero temperature IRLM | 28 | for the same hopping 
parameter. The solid curves are polynomial interpolations as 
guide to eye. 



t 


FIG. 3. Current in HF approximation as function of time 
for the IRLM for U' = 2to, t' = 0.5to, and several voltages 
close to the sharp transition at 14 ~ 1.9to/e in the I-V curve 
in Fig.[^ 


Ising variables, each of them representing the time evo¬ 
lution of a non-interacting system in a time-dependent 
potential. For not too long leads and a number of time 
slices M < 10 one can compute these sums numerically 
without introducing any further source of error. The 
discrete Hubbard-Stratonovich decoupling has previously 
been used in combination with an iterative summation 
scheme in order to calculate transport properties of the 
single impurity Anderson model uni. We will use the 
Hubbard-Stratonovich method later to benchmark the 
results of the Hartree-Fock approximation. 


RESULTS 

In the following we present numerical results for the 
model defined in Eq. ([^. We restrict ourselves to the 
cases Nq = 1 which corresponds to the IRLM, and 
Nc = 2 which we will refer to as two-site model, for 
brevity. In the IRLM U' is the only interaction parame¬ 
ter, whereas in the two-site model we set U' = 0 and vary 
the interaction U between the two atoms of the molecule. 

The interaction dependence of the linear conductance 
for a two-site model similar to ours was studied in |25] 
using DMRG. In contrast to our work, in that paper the 
interaction between molecule and leads was varied, while 
the interaction on the molecule was kept fixed. 

Unless stated otherwise, in all calculations the length 
of the leads is Afp = N-^ = 100, the inverse temperature 
is /3l = /3r = 20/to, and the overall chemical potential 
of the unbiased system is fj, = 0 which corresponds to 
half-filled band in the leads. The time-dependent cur¬ 
rent / is calculated as the mean value of left and right 
currents if they are different. The currents obtained from 
the plateau value reached in the time evolution are iden¬ 
tical to the steady state currents calculated from the self- 
consistent HF approach within a relative error typically 
of the order of 10“®, with the exception of the cases where 
hysteresis occurs in the self-consistent solution. 


Nc = 1: Interacting resonant level model 

In the case where the molecule consists of a single site 
the Hamiltonian 0 is identical to the IRLM with in¬ 
teraction parameter U'. The I-V characteristics of this 
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model for the special value U' = 2to is known analytically 
[ 20 ] . and is given in closed form by [26] 


e^V 
= Mi 


1 ? 1 

4’ 4’ 


5 7 
6 ’ 6 


(19) 

where 3 F 2 is a hypergeometric function, and eVcjto = 
The prefactor r « 3.2 is the lattice regu¬ 
larization of the corresponding field theory. In Fig. 2 
we compare the I-V characteristics obtained within the 
HF approximation with the exact analytical result of 


Eq. (19) for hopping parameters f = 0.3to and t' = OMq. 


For small voltages there is a linear relation between cur¬ 
rent and voltage due to the fact that at half filling the 
Fermi level is exactly at the resonance, and therefore 
the conductance is identical to the conductance quan¬ 
tum Go = e^//i of spinless fermions. Expanding Eq. (19) 
to leading order in the voltage yields 


( 20 ) 


/(v')=.G.r 



For comparison, the relative deviation A of the HF cur¬ 
rents from perfect conduction is shown in the inset of 
Fig. 2. A linear ht of the data on a double logarithmic 
scale yields A oc just like for the non-interacting 
system, in contrast to the nontrivial analytical result 
A oc E®. Perturbative calculations [57] for the IRLM 
away from the self-dual point U' = 2to indicate that 
there exists a nonvanishing contribution A cx E^ to the 
backscattered nonlinear conductance although with a re¬ 
duced weight compared to the noninteracting system. 

For larger voltages, the exact currents reach a maxi¬ 
mum and then decrease slowly but steadily with negative 
differential conductance. The HF currents, on the other 
hand, increase up to somewhat larger voltages. Then 
a sudden drop occurs followed by a range of negative 
differential conductance quite close to the analytical re¬ 
sult. For voltages beyond the band width, eE > 4<o, the 
HF current remains constant while the exact one contin¬ 
ues to decrease. The location of the jump is indepen¬ 
dent of the method how the HF currents are calculated. 
In particular there is no sign of hysteresis, i.e., solving 
the self-consistency equations for voltages increasing by 
small steps by taking the converged result of the previous 
voltage as starting point for the next iteration leads to 
the same curves as for incrementally decreasing voltages. 
The discontinuous behavior, which is obviously an arti¬ 
fact of the HF approximation, can also be observed in 
the time-evolution of the current displayed in Fig. 3 for 
t' = 0.5to and several voltages close to the jump. While 
for voltages below the transition (red curves) the sta¬ 
tionary state is reached quite fast and the plateau values 
increase with voltage, there is a quite strong reduction of 
the stationary current within a very small range of volt¬ 
ages, and it takes much longer for the system to reach 
the non-equilibrium steady state. 



FIG. 5. HF transmission coefficient of the two-site model at 
zero voltage for t' = O.Sto and U'/to = —1, 0,1,1.5 (from left 
to right). 


Figure 4 shows the HF transmission coefficient T{E) 
of the unbiased (zero voltage) IRLM which is closely re¬ 
lated to the local spectral function A{E). At half Filing 
A{E) oc Af{E)T{E), with the density of states 


The most striking feature is a pronounced broadening of 
the central peak with increasing interaction. The zero 
temperature spectral function of the IRLM has been cal¬ 
culated numerically using the Chebyshev expansion [28| , 
and the broadening of the central peak has been obtained 
by fitting the data to a Lorentzian, 

In the inset of Fig. 4 we compare our data with the broad¬ 
ening calculated from the exact spectral function. Up 
to U ~ 0.5to the HF results agree reasonably well with 
the exact ones but for larger interactions there is a quite 
strong discrepancy indicating the limitations of the HF 
approximation. 



Nc = 2: Two-site model 

The transport properties of the model ([^ in the case 
Nc = 2 are different from those of the IRLM in many 
respects due to the fact that now the Fermi level lies 
exactly between the two transmission resonance peaks of 
the non-interacting system. Therefore, the inclusion of 
interaction may not only broaden or shift these peaks but 
also modify the transmission at small energies and thus 
strongly influence the conductance. 

Figure 5 shows the HF transmission coefficient of the 
unbiased two-site model for several values of the interac¬ 
tion U. While for attractive interaction, U = —to, the 
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FIG. 6. Current voltage characteristics of the two-site model 
for t' — 0.5to, P = 8/to, and several values of U. The curves 
are HF data, the symbols are obtained from the Hubbard 
Stratonovich approach for shorter leads with A^l = Nb. = 20, 
as illustrated in Fig. 7. The size of the symbols is larger than 
the estimated error. The dotted curve is the HF current for 
U — 1.5to and /3 = 50/to. 


resonance peak is shifted to smaller energy values and 
slightly broadened with respect to the non-interacting 
one, the opposite is the case for repulsive interaction. 
Correspondingly, the transmission at zero energy is re¬ 
duced with increasing U which is expected to result in a 
decreased conductance. This expectation is confirmed 
in the I-V diagram (Fig. 6) where the linear conduc¬ 
tance is largest for attractive interaction U = —tg, and 
decreases with increasing values of U. Remarkably, at 
higher voltages there is a strong and for U > 1.5to even 
discontinuous increase of the current such that all curves 
nearly collapse in the large voltage regime. 

In order to check if the strong increase of the currents 
in the I-V characteristics is an artifact of the HF approx¬ 
imation, we have calculated the exact time-dependent 
currents using the Hubbard-Stratonovich approach de¬ 
scribed in section 2.4 for shorter leads (A^l = Nr = 20). 
The number of Trotter time slices M = 9 is chosen such 
that discretization errors in the data shown in the fig¬ 
ure are much smaller than the size of the symbols. To 
avoid finite size effects due to the discrete spectra of the 
leads, we have used a somewhat smaller inverse temper¬ 
ature, /3 = 8/to, for the data shown in Fig. 6 and in 
Fig. 7. To get an idea about the influence of finite tem¬ 
peratures we also show the HF current for U = 1.5to and 
/3 = 50/to as dotted curve in Fig. 6. The deviation from 
the /3 = 8/to result is negligible with the exception of 
the voltage region close to the jump. For small voltages, 
V < 2to/e, the HF currents agree reasonably well with 
the exact ones displayed as open symbols of the same 
color, whereas for larger voltages they are completely off. 
To elucidate this behavior, the time-dependent HF cur¬ 
rents are compared with the exact ones in Fig. 7 for two 



FIG. 7. Time-dependent currents for the two-site model with 
A'’l = A^’r = 20, P = 8/to, t' = 0.5to, and several values of 
the interaction strength. The symbols are the exact currents 
from the Hubbard Stratonovich approach for M = 9 time 
slices, the full curves are HF data. The arrows indicate the 
HF stationary state currents for the same parameters in the 
limit of infinitely long leads. Top panel: V = to/e, bottom 
panel: V = Ato/e. 

different voltages, V = to/e and V = Atg/e. For the 
smaller voltage (upper panel), the HF currents for re¬ 
pulsive interaction nearly coincides with the exact cur¬ 
rents and deviate only slightly for repulsive interaction, 
U = —tg. Note that the stationary currents for much 
longer leads indicated by arrows on the right axis are 
nearly identical to what one would obtain by averaging 
over the small oscillations observed in the HF data for 
Nl = 20 . 

In the case of large voltage (lower panel), the HF cur¬ 
rents resemble the exact ones only for short times during 
the transient phase of the time evolution but converge 
all to the same current plateau of the non-interacting 
system. The exact currents on the other hand approach 
different stationary values depending on the interaction 
strength U. 

It is therefore clear that the collapse of the different 
curves at large voltages in the I-V diagram of Fig. 6 is 
an artifact of the HF approach and indicates a fundamen- 
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FIG. 8. Hysteresis in the current voltage characteristics of 
the two-site model for U = l.Sto, t' = O.Sto, and P = SO/to- 
The blue (red) cur ve is obtained from the iterative solution 
of the HF equation (131 for adiabatically increasing (decreas¬ 
ing) voltage, as indicated by the arrows. The black curve 
corresponds to the stationary current taken from the time 
evolution. 


yield identical results with comparable numerical effort. 
However, the self-consistent approach sometimes allows 
for multiple solutions which leads to hysteretic behavior 
when the voltage is varied adiabatically. This ambigu¬ 
ity can be avoided using the stationary current obtained 
from the time-evolution approach. For a model of in¬ 
teracting spinless fermions, the HF data agree well with 
available exact results, with the exception of the large 
voltage regime of the two-site model where a spurious 
discontinuous transition is observed within the HF ap¬ 
proximation. It is straightforward to generalize the HF 
approach in many respects. In addition to the charge cur¬ 
rents also energy or heat currents can be calculated which 
is of particular interest when besides the voltage there 
is also a temperature gradient. Furthermore, dynamical 
properties, e.g., the response to a time-dependent gate 
voltage, can be studied without significant additional ef¬ 
fort. 
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The existence of multiple solutions in the I-V charac¬ 
teristics within the HF approximation and within the 
adiabatic local density approximation of density func¬ 
tional theory has recently been discussed for a model with 
Hubbard-type interaction [2^ . 

CONCLUSION 

The time-dependent HF approximation is a compu¬ 
tationally cheap and versatile approach to calculate the 
I-V characteristics of weakly correlated systems at finite 
temperatures. The time-evolution of the currents until a 
plateau value is reached as well as an iterative solution of 
the self-consistent HF equations for the stationary state 
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